DESeq2 differential analysis
M Jaritz
2025-10-14
#remotes::install_github("thomasp85/patchwork")
#remotes::install_github("tidyverse/ggplot2", ref = remotes::github_pull("5592"))
#library(optparse)
#library(readr) tidyverse
#library(stringr) tidyverse
# R 4.3.2
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
library(tidyverse)
library(GGally)
library(ggrepel)
library(DESeq2)
library(edgeR) # cpm
library(vsn) # meanSdplot
library(ReportingTools) # HTMLReport
library(ape)
library(ShortRead)
library(apeglm)
library(EDASeq) # plotRLE
library(mixOmics)
library(matrixStats)
library(RColorBrewer)
library(viridis)
library(patchwork)
library(ggpointdensity)
library(data.table)
library(DT)
library(htmlTable)
library(DGEobj.utils)
source("/users/jaritz/Gaidt/bioinf/software/scripts/deseq2.plot.R")
source("/users/jaritz/Gaidt/bioinf/software/scripts/multiple_matrix_DESeq2.functions.R") # mPCA
source("/users/jaritz/Gaidt/bioinf/software/scripts/analyse.functions.R") # write_extended_results_tableworkdate <- "2025Oct13"
workdir <- paste0("/groups/gaidt/user/markus.jaritz/projects/Markus/",workdate,"/")
sampleTable=paste0(workdir,"doc/samples.txt")
dir.create(paste0(workdir,"/results"),recursive = T,showWarnings = F)
setwd(paste0(workdir,"/results"))
getwd()
odir <- getwd()
params=list(
sampleTable='/users/jaritz/Gaidt/bioinf/software/danbing-tk/samples_experiment.txt',
conditions_colname='condition',
control_condition='exp1',
experiment_condition='exp2',
countTableAnnot='/scratch-cbe/users/markus.jaritz/Gaidt/Moritz/2025Oct13/results/featureCount_concat/featureCount_counts.annot.txt',
countTable='/scratch-cbe/users/markus.jaritz/Gaidt/Moritz/2025Oct13/results/featureCount_concat/featureCount_counts.txt',
odir=odir,
downsample_pairs=0,
downsample_ma=0,
session="http://localhost:60151/goto?"
)
params
p_ <- GGally::print_if_interactive
set.seed(42)p_ <- GGally::print_if_interactive
set.seed(42)
print(paste("sampleTable: ",params$sampleTable))[1] "sampleTable: samples_experiment.txt"print(paste("conditions_colname: ",params$conditions_colname))[1] "conditions_colname: condition"print(paste("control_condition: ",params$control_condition))[1] "control_condition: exp1"print(paste("experiment_condition: ",params$experiment_condition))[1] "experiment_condition: exp2"print(paste("countTable: ",params$countTable))[1] "countTable: /groups/gaidt/user/markus.jaritz/projects/Markus/2025Oct13/results/featureCount_concat/featureCount_counts.txt"print(paste("countTableAnnot: ",params$countTableAnnot))[1] "countTableAnnot: /groups/gaidt/user/markus.jaritz/projects/Markus/2025Oct13/results/featureCount_concat/featureCount_counts.annot.txt"print(paste("downsample_pairs: ",params$downsample_pairs))[1] "downsample_pairs: 1000"print(paste("downsample_ma: ",params$downsample_ma))[1] "downsample_ma: 0"# DEBUG
# params <- list(celltype="BLaER1")1 INPUT
#
# count table
#
t <- read_delim(params$countTable,delim="\t",col_names = T,skip = 0)
print("input countTable:")[1] "input countTable:"print(dim(t))[1] 22246 5t[1:3,1:3]# A tibble: 3 × 3
sites `MORC3-V5_DMSO_M3AID_REP1_featureCount` `MORC3-V5_DMSO_M3AID_REP2_featureCount`
<chr> <dbl> <dbl>
1 site_1 44 26
2 site_2 72 47
3 site_3 26 13t[t$sites=="site_1468",]# A tibble: 1 × 5
sites `MORC3-V5_DMSO_M3AID_REP1_featureCount` `MORC3-V5_DMSO_M3AID_REP2_featureCount` sgAAVS1_DMSO_R1_featureCount sgAAVS1_DMSO_R2_featureCount
<chr> <dbl> <dbl> <dbl> <dbl>
1 site_1468 29 14 10 22#
# sample meta info
#
m <- read_delim(params$sampleTable,delim="\t",col_names = T)
m# A tibble: 4 × 4
sample_id sample_name count_column condition
<dbl> <chr> <chr> <chr>
1 1 MORC3_V5_DMSO_M3AID_REP1_exp2 MORC3-V5_DMSO_M3AID_REP1_featureCount exp2
2 2 MORC3_V5_DMSO_M3AID_REP2_exp2 MORC3-V5_DMSO_M3AID_REP2_featureCount exp2
3 3 sgAAVS1_DMSO_R1_exp1 sgAAVS1_DMSO_R1_featureCount exp1
4 4 sgAAVS1_DMSO_R2_exp1 sgAAVS1_DMSO_R2_featureCount exp1 datatable(m)#
# annotation data
#
countTableAnnot <- read_delim(params$countTableAnnot,delim="\t",col_names = T)
colnames(countTableAnnot)[1] <- "SiteID"
countTableAnnot[1:3,]# A tibble: 3 × 6
SiteID Chr Start End Strand Length
<chr> <chr> <dbl> <dbl> <chr> <dbl>
1 site_1 chr1 826913 827797 + 885
2 site_2 chr1 904446 905303 + 858
3 site_3 chr1 906294 907109 + 816countTableAnnot %>% filter(SiteID=="site_1468")# A tibble: 1 × 6
SiteID Chr Start End Strand Length
<chr> <chr> <dbl> <dbl> <chr> <dbl>
1 site_1468 chr1 146472327 146472590 + 264# counts only
tc <- t[,-1]
print(dim(tc))[1] 22246 4# setup DESeqDataSet
tcm <- as.matrix(tc)
tcm[1:3,1:3] MORC3-V5_DMSO_M3AID_REP1_featureCount MORC3-V5_DMSO_M3AID_REP2_featureCount sgAAVS1_DMSO_R1_featureCount
[1,] 44 26 60
[2,] 72 47 88
[3,] 26 13 76counts <- matrix(tcm,
ncol=dim(tcm)[2],
dimnames=list(t$sites,colnames(tcm)))
print("Final deseq2 count table:")[1] "Final deseq2 count table:"print(dim(counts))[1] 22246 4counts[1:3,] MORC3-V5_DMSO_M3AID_REP1_featureCount MORC3-V5_DMSO_M3AID_REP2_featureCount sgAAVS1_DMSO_R1_featureCount sgAAVS1_DMSO_R2_featureCount
site_1 44 26 60 69
site_2 72 47 88 88
site_3 26 13 76 392 DESeq2
DDS <- DESeqDataSetFromMatrix(countData = counts,
colData = m,
design = ~condition)
head(colData(DDS))DataFrame with 4 rows and 4 columns
sample_id sample_name count_column condition
<numeric> <character> <character> <factor>
MORC3-V5_DMSO_M3AID_REP1_featureCount 1 MORC3_V5_DMSO_M3AID_.. MORC3-V5_DMSO_M3AID_.. exp2
MORC3-V5_DMSO_M3AID_REP2_featureCount 2 MORC3_V5_DMSO_M3AID_.. MORC3-V5_DMSO_M3AID_.. exp2
sgAAVS1_DMSO_R1_featureCount 3 sgAAVS1_DMSO_R1_exp1 sgAAVS1_DMSO_R1_feat.. exp1
sgAAVS1_DMSO_R2_featureCount 4 sgAAVS1_DMSO_R2_exp1 sgAAVS1_DMSO_R2_feat.. exp1rowData(DDS) <- countTableAnnot
head(rowData(DDS))DataFrame with 6 rows and 6 columns
SiteID Chr Start End Strand Length
<character> <character> <numeric> <numeric> <character> <numeric>
site_1 site_1 chr1 826913 827797 + 885
site_2 site_2 chr1 904446 905303 + 858
site_3 site_3 chr1 906294 907109 + 816
site_4 site_4 chr1 921152 921458 + 307
site_5 site_5 chr1 924657 924957 + 301
site_6 site_6 chr1 940863 941872 + 1010print(dim(DDS))[1] 22246 4rowData(DDS)[rowData(DDS)$SiteID=="site_1468",]DataFrame with 1 row and 6 columns
SiteID Chr Start End Strand Length
<character> <character> <numeric> <numeric> <character> <numeric>
site_1468 site_1468 chr1 146472327 146472590 + 2642.1 Relative fold changes and PCA
The Relative Log Expression (RLE) plot is a useful diagnostic plot to visualize the differences between the distributions of read counts across samples.
It shows the boxplots of the log-ratios of the site-level read counts of each sample to those of a reference sample (defined as the median across the samples). Ideally, the distributions should be centered around the zero line and as tight as possible. Clear deviations indicate the need for normalization and/or the presence of outlying samples.
The PCA plot represents a Principal Component Analysis (PCA) plot of the raw counts in object .
#
# count check
#
par(las=2,mar=c(10,2,2,2))
plotRLE(counts(DDS))
plotPCA(counts(DDS))
#
# count distribution
#
# all counts of all samples together
raw_count_cuttoff <- 0
as_tibble(counts(DDS)) %>%
pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
ggplot(aes(x=log2(counts+1))) +
geom_histogram(bins=100) +
geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
theme_light() +
ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))
raw_count_cuttoff <- 1
as_tibble(counts(DDS)) %>%
pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
ggplot(aes(x=log2(counts+1))) +
geom_histogram(bins=100) +
geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
theme_light() +
ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))
# by sample
as_tibble(counts(DDS)) %>%
pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
ggplot(aes(x=log2(counts+1))) +
geom_histogram(bins=100) +
geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
facet_wrap(~samples) +
theme_light() +
ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))
# rowsums
rs <- tibble(raw_count_sums_per_peak=rowSums(counts(DDS)))
rs %>% ggplot(aes(x=log2(raw_count_sums_per_peak+1))) +
geom_histogram(bins=100) +
theme_light() +
ggtitle("Raw count rowsums (over all replicates, per site)")
#
# filter a bit for very low numbers
#
print("low counts filtered out:")[1] "low counts filtered out:"smallestGroupSize <- 2
keep <- rowSums(counts(DDS) > (2*raw_count_cuttoff)) >= smallestGroupSize
print(table(keep))keep
FALSE TRUE
2 22244 DDS <- DDS[keep,]
dim(DDS)[1] 22244 4plotRLE(counts(DDS))
plotPCA(counts(DDS))
sections=c(10,40,800,2000,5000)
plot_counts(counts(DDS)+1,sections)[1] "plotting bargraph ..."
[1] "plotting overall count sum ditribution ..."
[1] "plotting overall count sum ditribution, by time and replicate ..."
[1] "final plot ..."
2.2 RLOG transform (VST)
#
# get rlogs
#
if(file.exists(paste0(params$odir,"/vst.RData"))) {
print(paste0("loading ",params$odir,"/vst.RData"))
load(paste0(params$odir,"/vst.RData"))
} else {
print("calculating VST ...")
vst <- rlog(DDS, blind=T)
print(dim(DDS))
print(dim(vst))
save(vst,file=paste0(params$odir,"/vst.RData"))
}[1] "calculating VST ..."
[1] 22244 4
[1] 22244 42.2.1 PCA
We investigate the relationship between mean and variance, and its effect on the PCA.
time.var.png
mPCA(m=assay(vst),
colgroups=m$condition,
collegend="Condition",
pchgroups=as.character(m$sample_id),
pchlegend="Replicate",
#annot=rowData(vst)[,1],
title="vst.blind",
ofile=paste0(params$odir,"/",params$condition),
varCutoff = 0.9,
show_biplot=F)[1] "mPCA ..."
[1] "make matrix ..."
[1] 22244 4
[1] "pchgroups set to:"
[1] "1" "2" "3" "4"
[1] "pchgroups set to:"
1 2 3 4
1 2 3 4
[1] "calculate variance ..."
[1] "plotting variances ..."
[1] "violin ..."
[1] "density ..."
[1] "plotting variances versus means ..."
[1] "density ..."
[1] "scatter ..."
[1] "patchwork ..."[1] "selected features (variances only):"
[1] 2225
[1] "top selected features (variances only):"
NULL
[1] "taking default labels ...."
[1] "labels:"
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount" "MORC3-V5_DMSO_M3AID_REP2_featureCount" "sgAAVS1_DMSO_R1_featureCount" "sgAAVS1_DMSO_R2_featureCount"
[1] "pca ..."
[1] "plotIndiv ..."
[1] "colgroups_factor:"
[1] "exp2" "exp2" "exp1" "exp1"
[1] "pchgroup:"
1 2 3 4
1 2 3 4
[1] "pchlevels:"
[1] "1" "2" "3" "4"
[1] "selected features (variance and means):"
[1] "0.0703069724604185 5.35432301347482"
[1] 1214
[1] "top selected features (variance and means):"
NULL
[1] "selected features (no filtering, only variance == 0 filtered out):"
[1] 22244
[1] "top selected features (no filtering, only variance == 0 filtered out):"
NULL
[1] "selected features (variance>0):"
[1] 22244
2.3 Pairwise scatters RAW & VST
#
# plot for each celltype
#
# subsample, for quicker drawing
if(params$downsample_pairs>0) {
si <- sample(c(1:dim(DDS)[1]),params$downsample_pairs)
print(head(si))
DDS.sub <- DDS[si,]
} else {
DDS.sub <- DDS
si <- c(1:dim(DDS)[1])
}[1] 18753 21657 9290 1252 15506 8826# plot raw counts
p <- ggpairs(data.frame(log2(counts(DDS.sub)+1)),
lower=list(continuous=my_fn)) +
ggtitle(paste0("Raw counts subsample: ",params$downsample_pairs)) +
theme_bw()
#p_(p)
print(p)
ggsave(p,filename=paste0(params$odir,"/raw.ss_",params$downsample_pairs,".jpeg"))
# plot vst counts
vst.sub <- vst[si,]
p <- ggpairs(data.frame(assay(vst.sub)),
lower=list(continuous=my_fn)) +
ggtitle(paste0("RLOG transformed counts, subsample: ",params$downsample_pairs)) +
theme_bw()
#p_(p)
print(p)
ggsave(p,filename=paste0(params$odir,"/vst.ss_",params$downsample_pairs,".jpeg")) 2.4 Pairwise scatters TPM
#
# plot for each celltype
#
# subsample, for quicker drawing
if(params$downsample_pairs>0) {
si <- sample(c(1:dim(DDS)[1]),params$downsample_pairs)
print(head(si))
DDS.sub <- DDS[si,]
} else {
DDS.sub <- DDS
si <- c(1:dim(DDS)[1])
}[1] 7580 6620 4288 16478 5322 14895# plot TPMs
rd <- rowData(DDS.sub)
sitelengths <- rd$End-rd$Start
tpm.counts <- convertCounts(countsMatrix = counts(DDS.sub),unit="TPM",geneLength=sitelengths)
p <- ggpairs(data.frame(log2(tpm.counts+1)),
lower=list(continuous=my_fn)) +
ggtitle(paste0("TPMs, subsample: ",params$downsample_pairs)) +
theme_bw()
#p_(p)
print(p)
ggsave(p,filename=paste0(params$odir,"/tpm.ss_",params$downsample_pairs,".jpeg")) 2.5 Dispersion and normalized counts
plotDispEsts plots the per-gene dispersion estimates together with the fitted mean-dispersion relationship.
#
# ---- DESeq2 model ----
#
#if(file.exists(paste0(params$odir,"/deobj.",params$celltype,",RData"))) {
# print(paste0("loading ",params$odir,"/deobj.",params$celltype,".RData"))
# load(paste0(params$odir,"/deobj.",params$celltype,".RData"))
#} else {
print("calculating DESeq2 object ...")[1] "calculating DESeq2 object ..." deobj <- DESeq(DDS, test='Wald')
save(deobj,file=paste0(params$odir,"/deobj.RData"))
#}
print("DESeq2 sizeFactors:")[1] "DESeq2 sizeFactors:"p <- tibble("names"=names(sizeFactors(deobj)),
"sf"=sizeFactors(deobj)) %>%
ggplot(aes(x=names,y=sf)) +
geom_col() +
theme_light() +
theme(axis.text.x=element_text(angle=90)) +
ggtitle("Size factors")
print(p)
ggsave(p,filename=paste0(params$odir,"/sizefactors.jpeg"))
#
print("Dispersion estimates:")[1] "Dispersion estimates:"plotDispEsts(deobj)
res <- results(deobj)2.6 Relative fold change and PCA - normalized
par(las=2,mar=c(10,2,2,2))
plotRLE(counts(deobj,normalized=T))
plotPCA(counts(deobj,normalized=T))
2.7 Pairwise scatters NORM
#
# ---- DESeq2 plots ----
#
# ---- norm counts ----
cts.norm <- counts(deobj,normalized=T)
# same si as above
if(params$downsample_pairs>0) {
cts.norm.sub <- cts.norm[si,]
} else {
cts.norm.sub <- cts.norm
}
# summary(is.na(cts.norm)) # there should be no NA
p <- ggpairs(log2(cts.norm.sub+1),
lower=list(continuous=my_fn)) +
ggtitle(paste0("Normalized counts, subsample=",params$downsample_pairs)) +
theme_bw()
#p_(p)
print(p)
ggsave(p,filename=paste0(params$odir,"/norm.ss_",params$downsample_pairs,".jpeg")) 2.8 MA and volcano plots
if(params$downsample_ma>0) {
si <- sample(c(1:dim(deobj)[1]),params$downsample_ma)
deobj.sub <- deobj[si,]
vst.sub <- vst[si,]
} else {
si <- c(1:dim(deobj)[1])
deobj.sub <- deobj
vst.sub <- vst
}
res <- results(deobj.sub,contrast = c("condition",params$experiment_condition,params$control_condition))
res.df <- results(deobj.sub,contrast = c("condition",params$experiment_condition,params$control_condition),tidy=T)
dim(res.df)[1] 22244 7res.df[1:3,1:3] row baseMean log2FoldChange
1 site_1 46.75372 0.04882633
2 site_2 71.95839 0.37476977
3 site_3 33.45225 -0.62945101res.df[1:3,] row baseMean log2FoldChange lfcSE stat pvalue padj
1 site_1 46.75372 0.04882633 0.5478855 0.08911775 0.9289883 0.9994137
2 site_2 71.95839 0.37476977 0.4462411 0.83983702 0.4009998 0.9994137
3 site_3 33.45225 -0.62945101 0.6765031 -0.93044807 0.3521391 NA#data.frame(rowData(deobj.sub))
row.df <- data.frame(rowData(deobj.sub)[1:19])
row.df <- cbind(row=rownames(row.df),row.df)
row.df[1:3,1:3] row SiteID Chr
site_1 site_1 site_1 chr1
site_2 site_2 site_2 chr1
site_3 site_3 site_3 chr1res.df.a <- merge(x=row.df,y=res.df,by="row",sort=T,all=T)
res.df.a[1:3,] row SiteID Chr Start End Strand Length baseMean.x baseVar allZero dispGeneEst dispGeneIter dispFit dispersion dispIter dispOutlier dispMAP Intercept condition_exp2_vs_exp1 SE_Intercept baseMean.y
1 site_1 site_1 chr1 826913 827797 + 885 46.75372 39.33170 FALSE 0.00000001 2 0.1462236 0.1218074 6 FALSE 0.1218074 5.524677 0.04882633 0.3779915 46.75372
2 site_10 site_10 chr1 999572 1000802 + 1231 50.08206 454.39487 FALSE 0.26045883 3 0.1373907 0.1567706 9 FALSE 0.1567706 5.755564 -0.21399560 0.4205713 50.08206
3 site_100 site_100 chr1 2953145 2953475 + 331 10.95914 80.33834 FALSE 1.81124086 7 0.5803305 0.6815258 9 FALSE 0.6815258 3.980974 -1.32309113 0.8696165 10.95914
log2FoldChange lfcSE stat pvalue padj
1 0.04882633 0.5478855 0.08911775 0.9289883 0.9994137
2 -0.21399560 0.6093170 -0.35120570 0.7254340 0.9994137
3 -1.32309113 1.3021971 -1.01604523 0.3096078 NAdim(res.df.a)[1] 22244 26dim(vst.sub)[1] 22244 4dim(counts(deobj.sub))[1] 22244 4desc <- "test"
write_extended_results_table(deobj=deobj.sub,
vst=vst.sub,
res=res.df.a,
ofile=paste0(params$odir,"/",desc,".",params$downsample_ma,".txt"),
genelength = res.df.a$End-res.df.a$Start)[1] 22244
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244 4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244 4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244 4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"#
# MA plots
#
reslog2 fold change (MLE): condition exp2 vs exp1
Wald test p-value: condition exp2 vs exp1
DataFrame with 22244 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
site_1 46.7537 0.04882633 0.547886 0.08911775 0.928988 0.999414
site_2 71.9584 0.37476977 0.446241 0.83983702 0.401000 0.999414
site_3 33.4523 -0.62945101 0.676503 -0.93044807 0.352139 NA
site_4 21.0182 -0.00537712 0.801507 -0.00670876 0.994647 NA
site_5 16.1062 -0.67071868 0.949427 -0.70644553 0.479911 NA
... ... ... ... ... ... ...
site_22242 19.6523 -0.331930 0.846293 -0.392217 0.694898 NA
site_22243 17.4983 -1.117983 0.888502 -1.258279 0.208291 NA
site_22244 49.2172 -0.271975 0.615693 -0.441737 0.658679 0.999414
site_22245 19.8220 -0.191921 0.826885 -0.232101 0.816460 NA
site_22246 30.5679 0.191471 0.702078 0.272720 0.785069 NAp <- my_ma(padj_co=0.01,deobj=res)
print(p)NULLp1 <- my_ma(padj_co=0.01,deobj=NULL,df=data.frame(res),baseMean_co=0, color=T,
filter_df_na = T, name=desc,
outfile=paste0(params$odir,"/ma.",params$downsample_ma,".jpeg"))
print(p1)
p2 <- my_volcano(df=data.frame(res),baseMean_co = 0,filter_df_na = T,
color=T, name="test",
outfile=paste0(params$odir,"/vc.",params$downsample_ma,".jpeg"))
print(p2)
# top results table with links
session="http://localhost:60151/goto?"
pad=10000
html.df <- coords_to_igv_links(df = res.df.a %>% arrange(pvalue) %>% slice_head(n=100),
chr="Chr",start="Start",end="End",session=session,pad=pad) # params$session)
html.table <- html.df %>%
addHtmlTableStyle(align = "r") %>%
htmlTable(caption=paste("test","top 100 entries, sorted by pvalue, pad=",pad),
useViewer=FALSE,rnames=FALSE)
write(html.table,file = paste0(params$odir,"/test.",params$downsample_ma,".top100.html"))3 Session info
sessionInfo()R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.11.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Vienna
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] DGEobj.utils_1.0.6 htmlTable_2.4.2 DT_0.31 data.table_1.15.4 ggpointdensity_0.1.0 patchwork_1.2.0 viridis_0.6.5
[8] viridisLite_0.4.2 RColorBrewer_1.1-3 mixOmics_6.26.0 lattice_0.22-6 MASS_7.3-60.0.1 EDASeq_2.36.0 apeglm_1.24.0
[15] ShortRead_1.60.0 GenomicAlignments_1.38.2 Rsamtools_2.18.0 Biostrings_2.70.3 XVector_0.42.0 BiocParallel_1.36.0 ape_5.8
[22] ReportingTools_2.42.3 knitr_1.45 vsn_3.70.0 edgeR_4.0.16 limma_3.58.1 DESeq2_1.42.1 SummarizedExperiment_1.32.0
[29] Biobase_2.62.0 MatrixGenerics_1.14.0 matrixStats_1.3.0 GenomicRanges_1.54.1 GenomeInfoDb_1.38.8 IRanges_2.36.0 S4Vectors_0.40.2
[36] BiocGenerics_0.48.1 ggrepel_0.9.5 GGally_2.2.1 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.4
[43] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] ProtGenerics_1.34.0 bitops_1.0-7 httr_1.4.7 Rgraphviz_2.46.0 numDeriv_2016.8-1.1 tools_4.3.2 backports_1.4.1 utf8_1.2.4
[9] R6_2.5.1 mgcv_1.9-1 lazyeval_0.2.2 withr_3.0.0 prettyunits_1.2.0 gridExtra_2.3 preprocessCore_1.64.0 textshaping_0.3.7
[17] cli_3.6.2 labeling_0.4.3 ggbio_1.50.0 sass_0.4.7 mvtnorm_1.2-4 genefilter_1.84.0 systemfonts_1.0.5 foreign_0.8-86
[25] R.utils_2.12.3 AnnotationForge_1.44.0 dichromat_2.0-0.1 BSgenome_1.70.2 bbmle_1.0.25.1 rstudioapi_0.15.0 RSQLite_2.3.6 generics_0.1.3
[33] GOstats_2.68.0 BiocIO_1.12.0 hwriter_1.3.2.1 crosstalk_1.2.1 vroom_1.6.5 GO.db_3.18.0 Matrix_1.6-5 interp_1.1-6
[41] fansi_1.0.5 abind_1.4-5 R.methodsS3_1.8.2 lifecycle_1.0.4 yaml_2.3.7 SparseArray_1.2.4 BiocFileCache_2.10.2 grid_4.3.2
[49] blob_1.2.4 crayon_1.5.2 bdsmatrix_1.3-7 GenomicFeatures_1.54.4 annotate_1.80.0 KEGGREST_1.42.0 pillar_1.9.0 rjson_0.2.21
[57] corpcor_1.6.10 codetools_0.2-20 glue_1.7.0 vctrs_0.6.5 png_0.1-8 gtable_0.3.5 assertthat_0.2.1 emdbook_1.3.13
[65] cachem_1.0.8 aroma.light_3.32.0 xfun_0.41 S4Arrays_1.2.1 coda_0.19-4.1 survival_3.6-4 statmod_1.5.0 ellipsis_0.3.2
[73] nlme_3.1-164 Category_2.68.0 bit64_4.0.5 progress_1.2.3 filelock_1.0.3 bslib_0.5.1 affyio_1.72.0 rpart_4.1.23
[81] colorspace_2.1-0 DBI_1.2.2 Hmisc_5.1-2 nnet_7.3-19 tidyselect_1.2.1 bit_4.0.5 compiler_4.3.2 curl_5.1.0
[89] graph_1.80.0 xml2_1.3.5 DelayedArray_0.28.0 rtracklayer_1.62.0 checkmate_2.3.1 scales_1.3.0 hexbin_1.28.3 affy_1.80.0
[97] RBGL_1.78.0 rappdirs_0.3.3 digest_0.6.33 rmarkdown_2.25 htmltools_0.5.7 pkgconfig_2.0.3 jpeg_0.1-10 base64enc_0.1-3
[105] highr_0.10 dbplyr_2.5.0 fastmap_1.1.1 ensembldb_2.26.0 rlang_1.1.3 htmlwidgets_1.6.2 PFAM.db_3.18.0 farver_2.1.1
[113] jquerylib_0.1.4 jsonlite_1.8.7 R.oo_1.26.0 VariantAnnotation_1.48.1 RCurl_1.98-1.14 magrittr_2.0.3 Formula_1.2-5 GenomeInfoDbData_1.2.11
[121] munsell_0.5.1 Rcpp_1.0.11 stringi_1.7.12 zlibbioc_1.48.2 plyr_1.8.9 ggstats_0.6.0 parallel_4.3.2 deldir_2.0-4
[129] splines_4.3.2 hms_1.1.3 locfit_1.5-9.9 igraph_2.0.3 reshape2_1.4.4 biomaRt_2.58.2 DGEobj_1.1.2 XML_3.99-0.16.1
[137] evaluate_0.23 latticeExtra_0.6-30 biovizBase_1.50.0 BiocManager_1.30.22 tzdb_0.4.0 xtable_1.8-4 restfulr_0.0.15 AnnotationFilter_1.26.0
[145] RSpectra_0.16-1 ragg_1.2.6 rARPACK_0.11-0 OrganismDbi_1.44.0 memoise_2.0.1 AnnotationDbi_1.64.1 ellipse_0.5.0 cluster_2.1.6
[153] timechange_0.3.0 GSEABase_1.64.0